Antenna reader example

from radiocalibrationtoolkit import *
from healpy.newvisufunc import projview
import plotly.graph_objects as go

# This ensures Plotly output works in multiple places:
# plotly_mimetype: VS Code notebook UI
# notebook: "Jupyter: Export to HTML" command in VS Code
# See https://plotly.com/python/renderers/#multiple-renderers
import plotly.io as pio
pio.renderers.default = "plotly_mimetype+notebook"

[INFO] LFmap: Import successful.

layout_settings = dict(
    xaxis=dict(title="azimuth [°]", tickprefix="<b>", ticksuffix="</b>", dtick=30),
    yaxis=dict(
        title="zenith angle [°]",
        tickprefix="<b>",
        ticksuffix="</b>",
        dtick=10,
    ),
    coloraxis=dict(
        colorbar=dict(
            title=dict(
                text="VEL",
                side="right",
            ),
            tickprefix="<b>",
            ticksuffix="</b>",
        ),
    ),
    font=dict(
        size=15,
        color="black",
    ),
)

colormap = 'jet'

def get_plotly_cbar_label(quantity):
    if quantity == 'absolute':
        return "<|H|>"
    elif "Phi_amp" in quantity:
        return "|H<sub>Φ</sub>|"
    elif "Theta_amp" in quantity:
        return "|H<sub>ϴ</sub>|"

# create antenna instance
antenna_inst = AntennaPattern("./antenna_setup_files/SALLA_EW.xml")

[INFO] Your keys are: ['EAHTheta_amp', 'EAHTheta_phase', 'EAHPhi_amp', 'EAHPhi_phase', 'absolute'] Use them as the 'quantity'

# XML values in dictionary
# antenna_inst.get_raw()

# get antenna gain for unpolarized emission
df = antenna_inst.get(frequency=45, quantity="absolute")

# plot
fig = px.imshow(df.T, width=600, aspect="cube", origin='lower', color_continuous_scale=colormap)
fig.update_layout(**layout_settings)
fig.update_layout(
    coloraxis=dict(colorbar=dict(title=dict(text=get_plotly_cbar_label("absolute"))))
)
fig.show()

# interpolate
df = antenna_inst.get(frequency=45, quantity='absolute', interp_phi=np.linspace(0, 360, 200), interp_theta=np.linspace(0,90, 100))

# plot
fig = px.imshow(df.T, width=600, aspect="cube", origin='lower', color_continuous_scale=colormap)
fig.update_layout(**layout_settings)
fig.update_layout(
    coloraxis=dict(colorbar=dict(title=dict(text=get_plotly_cbar_label("absolute"))))
)
fig.show()

# antenna gain as a function of frequency, azimuth and zenith in the required polarization as volumetric dataset
volumetric_dataset_df = antenna_inst.get_volumetric_dataset(
    quantity="EAHTheta_amp", frequencies=np.arange(30, 80, 1)
)
volumetric_dataset_df

	theta	0.0	1.5	3.0	4.5	6.0	7.5	9.0	10.5	12.0	13.5	...	76.5	78.0	79.5	81.0	82.5	84.0	85.5	87.0	88.5	90.0
freq	phi
30	0.0	1.60419	1.60520	1.60644	1.60790	1.60954	1.61132	1.61319	1.61509	1.61695	1.61869	...	1.04430	0.99178	0.92724	0.84915	0.75593	0.64597	0.51750	0.36859	0.19698	0.0
	3.0	1.60425	1.60522	1.60643	1.60786	1.60948	1.61124	1.61309	1.61497	1.61682	1.61855	...	1.04482	0.99226	0.92768	0.84954	0.75627	0.64625	0.51773	0.36874	0.19706	0.0
	6.0	1.59993	1.60086	1.60205	1.60345	1.60504	1.60678	1.60861	1.61048	1.61231	1.61402	...	1.04277	0.99030	0.92584	0.84784	0.75475	0.64495	0.51667	0.36799	0.19666	0.0
	9.0	1.59122	1.59213	1.59329	1.59467	1.59623	1.59795	1.59977	1.60162	1.60344	1.60514	...	1.03815	0.98591	0.92171	0.84405	0.75137	0.64205	0.51435	0.36634	0.19578	0.0
	12.0	1.57817	1.57905	1.58019	1.58154	1.58309	1.58478	1.58658	1.58842	1.59022	1.59191	...	1.03097	0.97908	0.91531	0.83818	0.74614	0.63758	0.51076	0.36378	0.19441	0.0
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
79	348.0	1.45920	1.46010	1.45994	1.45871	1.45636	1.45283	1.44810	1.44215	1.43499	1.42669	...	1.49450	1.53392	1.54726	1.51992	1.43914	1.29492	1.08041	0.79232	0.43093	0.0
	351.0	1.47382	1.47481	1.47469	1.47346	1.47107	1.46746	1.46263	1.45653	1.44920	1.44068	...	1.51165	1.55085	1.56376	1.53564	1.45367	1.30774	1.09096	0.79998	0.43507	0.0
	354.0	1.48440	1.48547	1.48540	1.48419	1.48178	1.47814	1.47323	1.46703	1.45956	1.45090	...	1.52568	1.56459	1.57702	1.54819	1.46519	1.31785	1.09924	0.80598	0.43831	0.0
	357.0	1.49092	1.49206	1.49205	1.49087	1.48847	1.48483	1.47988	1.47363	1.46609	1.45734	...	1.53642	1.57494	1.58688	1.55740	1.47356	1.32515	1.10518	0.81025	0.44062	0.0
	360.0	1.49334	1.49456	1.49462	1.49349	1.49113	1.48750	1.48257	1.47631	1.46875	1.45998	...	1.54367	1.58174	1.59317	1.56314	1.47866	1.32950	1.10867	0.81275	0.44195	0.0

6050 rows × 61 columns

# zenith slice
df = volumetric_dataset_df.unstack().stack(level=0).loc[:, 66, :]
# plot
fig = px.imshow(df, width=600, aspect="cube", origin='lower', color_continuous_scale=colormap)
fig.update_layout(**layout_settings)
fig.update_layout(
    title="Slice at 66° zenith angle",
    xaxis=dict(title="azimuth [°]"),
    yaxis=dict(title="frequency [MHz]"),
    coloraxis=dict(
        colorbar=dict(title=dict(text=get_plotly_cbar_label("EAHTheta_amp")))
    ),
)

fig.show()

# azimuth slice
df = volumetric_dataset_df.loc[:, 180, :]
# plot
fig = px.imshow(df, width=600, aspect="cube", origin='lower', color_continuous_scale=colormap)
fig.update_layout(**layout_settings)
fig.update_layout(
    title="Slice at 180° azimuth",
    xaxis=dict(title="zenith angle [°]", dtick=10),
    yaxis=dict(title="frequency [MHz]"),
    coloraxis=dict(
        colorbar=dict(title=dict(text=get_plotly_cbar_label("EAHTheta_amp")))
    ),
)
fig.show()

## show antenna pattern as volumetric data with slices in all axis
quantity = "EAHTheta_amp"
volumetric_dataset_df = antenna_inst.get_volumetric_dataset(
    quantity=quantity, frequencies=np.arange(30, 80, 3)
)

volume_data = volumetric_dataset_df.to_numpy().reshape(
    volumetric_dataset_df.index.get_level_values("freq").nunique(),
    -1,
    volumetric_dataset_df.shape[1],
)
PHI, FREQ, THETA = np.meshgrid(
    volumetric_dataset_df.index.get_level_values("phi").unique().values,
    volumetric_dataset_df.index.get_level_values("freq").unique().values,
    volumetric_dataset_df.columns.values,
)


# some settings
font = "Arial Black"
colorscale = "jet"

scene = dict(
    xaxis=dict(title="<b>Zenith angle [°]</b>", color="black", dtick="30"),
    yaxis=dict(title="<b>Azimuth [°]</b>", color="black", dtick="90"),
    zaxis=dict(
        title="<b>frequency [MHz]</b>",
        tickvals=[40, 50, 60, 70, 80],
    ),
    aspectratio=dict(x=0.75, y=0.75, z=0.75),  # Adjust the aspect ratio as needed
)


def quantity_label2colorbar_dict(quantity_label):
    return dict(
        len=0.75,
        title=dict(
            text="<b>" + quantity_label + " [m]</b>",
            side="right",
        ),
        tickprefix="<b>",
        ticksuffix="</b>",
    )


def create_volume_trace(
    opacity=1,
    surface_fill=0.001,
    surface_count=1,
    opacityscale="uniform",
    caps=dict(x_show=False, y_show=False, z_show=False),
):
    """
    Create a go.Volume trace based on the given parameters.

    Parameters
    ----------
    opacity : float, optional
        The opacity of the volume. Default is 1.
    surface_fill : float, optional
        The surface fill value. Default is 0.001.
    surface_count : int, optional
        The surface count value. Default is 1.
    opacityscale : str, optional
        The opacity scale. Default is "uniform".
    caps : dict, optional
        The caps configuration. Default is dict(x_show=False, y_show=False, z_show=False).

    Returns
    -------
    go.Volume
        A go.Volume trace object.
    """
    return go.Volume(
        y=PHI.flatten(),
        x=THETA.flatten(),
        z=FREQ.flatten(),
        value=np.around(volume_data.flatten(), 3),
        cmin=0,
        cmax=2,
        opacity=opacity,
        # isomax=1.3,
        surface_fill=surface_fill,
        surface_count=surface_count,
        colorscale="jet",
        opacityscale=opacityscale,
        caps=caps,
        colorbar=quantity_label2colorbar_dict(get_plotly_cbar_label(quantity)),
    )


# # create volume data
# PHI, FREQ, THETA, volume_data, quantity_label = create_volume(
#     antenna_inst, quantity="EAHTheta_amp"
# )

# define some slices
slice_types = {
    "no_slices": [],  # Neutral entry for no slices
    "slices_z": [45, 60],
    "slices_x": [30, 65],
    "slices_y": [0, 270],
}

# make plots
for slice_type, locations in slice_types.items():
    if slice_type == "no_slices":
        fig = go.Figure()
        fig.add_trace(
            create_volume_trace(
                opacity=0.8,
                surface_fill=1,
                opacityscale="max",
                surface_count=40,
                caps=dict(x_show=True, y_show=True, z_show=True),
            )
        )
    else:
        fig = go.Figure()
        slice_trace = create_volume_trace()
        slice_trace.update(**{slice_type: dict(show=True, locations=locations)})
        fig.add_trace(slice_trace)

    fig.update_layout(
        scene=scene,
        margin=dict(r=50, b=10, l=10, t=10),
        height=600,
        autosize=False,
        font=dict(family=font, size=18, color="black"),
    )
    fig.show()

## show antenna pattern as volumetric data with slices in all axis
# version of plots with sliders


def show_slices(x, y, z, volume_data, slice_type="x"):
    """
    Show slices of volumetric data using the given parameters.

    Parameters
    ----------
    x : np.ndarray
        The x-coordinate values.
    y : np.ndarray
        The y-coordinate values.
    z : np.ndarray
        The z-coordinate values.
    volume_data : np.ndarray
        The volumetric data.
    slice_type : str, optional
        The type of slice. Valid values are 'x', 'y', or 'z'. Default is 'x'.

    Returns
    -------
    None
        This function does not return anything. It displays the plot.

    """
    x_delta = np.diff(x)[0]
    y_delta = np.diff(y)[0]
    z_delta = np.diff(z)[0]

    cmin = 0
    cmax = 2
    colorscale = "jet"

    if slice_type == "x":
        slider_label = x
        fig = go.Figure(
            frames=[
                go.Frame(
                    data=go.Surface(
                        x=np.ones(x.size) * x[k],
                        y=y,
                        z=(np.ones((z.size, y.size)) * z[:, np.newaxis]).T,
                        surfacecolor=(volume_data[:, :, k]).T,
                        cmin=cmin,
                        cmax=cmax,
                    ),
                    name=str(k),
                )
                for k in range(x.size)
            ]
        )

        # default fig
        fig.add_trace(
            go.Surface(
                x=np.ones(x.size) * x[0],
                y=y,
                z=(np.ones((z.size, y.size)) * z[:, np.newaxis]).T,
                surfacecolor=(volume_data[:, :, 0]).T,
                colorscale=colorscale,
                cmin=cmin,
                cmax=cmax,
                colorbar=quantity_label2colorbar_dict(quantity_label),
            )
        )

    elif slice_type == "y":
        slider_label = y
        fig = go.Figure(
            frames=[
                go.Frame(
                    data=go.Surface(
                        x=x,
                        y=np.ones(y.size) * y[k],
                        z=np.ones((z.size, x.size)) * z[:, np.newaxis],
                        surfacecolor=(volume_data[:, y.size - 1 - k, :]),
                        cmin=cmin,
                        cmax=cmax,
                    ),
                    name=str(k),
                )
                for k in range(y.size)
            ]
        )

        # default fig
        fig.add_trace(
            go.Surface(
                x=x,
                y=np.ones(y.size) * y[0],
                z=np.ones((z.size, x.size)) * z[:, np.newaxis],
                surfacecolor=(volume_data[:, 0, :]),
                colorscale=colorscale,
                cmin=cmin,
                cmax=cmax,
                colorbar=quantity_label2colorbar_dict(quantity_label),
            )
        )
    elif slice_type == "z":
        slider_label = z

        fig = go.Figure(
            frames=[
                go.Frame(
                    data=go.Surface(
                        x=x,
                        y=y,
                        z=z[k] * np.ones((y.size, x.size)),
                        surfacecolor=np.flipud(volume_data[z.size - 1 - k]),
                        cmin=cmin,
                        cmax=cmax,
                    ),
                    name=str(k),
                )
                for k in range(z.size)
            ]
        )

        # default fig
        fig.add_trace(
            go.Surface(
                x=x,
                y=y,
                z=z[0] * np.ones((y.size, x.size)),
                surfacecolor=np.flipud(volume_data[z.size - 1]),
                colorscale=colorscale,
                cmin=cmin,
                cmax=cmax,
                colorbar=quantity_label2colorbar_dict(quantity_label),
            )
        )

    def frame_args(duration):
        return {
            "frame": {"duration": duration},
            "mode": "immediate",
            "fromcurrent": True,
            "transition": {"duration": duration, "easing": "linear"},
        }

    sliders = [
        {
            "pad": {"b": 10, "t": 60},
            "len": 0.9,
            "x": 0.1,
            "y": 0,
            "steps": [
                {
                    "args": [[f.name], frame_args(0)],
                    "label": str(slider_label[k]),
                    "method": "animate",
                }
                for k, f in enumerate(fig.frames)
            ],
        }
    ]

    # Layout
    fig.update_layout(
        title="<br>Slices in volumetric data",
        width=800,
        height=800,
        scene=dict(
            yaxis=dict(
                range=[y[0] - y_delta, y[-1] + y_delta], autorange=False
            ),  # Set the y-axis range from -5 to +5
            zaxis=dict(range=[z[0] - z_delta, z[-1] + z_delta], autorange=False),
            xaxis=dict(range=[x[0] - x_delta, x[-1] + x_delta], autorange=False),
            aspectratio=dict(x=1, y=1, z=1),
        ),
        updatemenus=[
            {
                "buttons": [
                    {
                        "args": [None, frame_args(50)],
                        "label": "&#9654;",  # play symbol
                        "method": "animate",
                    },
                    {
                        "args": [[None], frame_args(0)],
                        "label": "&#9724;",  # pause symbol
                        "method": "animate",
                    },
                ],
                "direction": "left",
                "pad": {"r": 10, "t": 70},
                "type": "buttons",
                "x": 0.1,
                "y": 0,
            }
        ],
        sliders=sliders,
    )

    fig.update_layout(
        scene=scene,
        margin=dict(r=50, b=10, l=10, t=10),
        font=dict(family=font, size=18, color="black"),
    )

    fig.show()


# create arrays
quantity = "EAHTheta_amp"
quantity_label = get_plotly_cbar_label(quantity)

volumetric_dataset_df = antenna_inst.get_volumetric_dataset(
    quantity=quantity, frequencies=np.arange(30, 80, 3)
)

volume_data = volumetric_dataset_df.to_numpy().reshape(
    volumetric_dataset_df.index.get_level_values("freq").nunique(),
    -1,
    volumetric_dataset_df.shape[1],
)
x = volumetric_dataset_df.columns.values
y = volumetric_dataset_df.index.get_level_values("phi").unique().values
z = volumetric_dataset_df.index.get_level_values("freq").unique().values

# create slices with sliders
show_slices(x, y, z, volume_data, slice_type="x")
show_slices(x, y, z, volume_data, slice_type="y")
show_slices(x, y, z, volume_data, slice_type="z")

# convert to healpy format
update_antenna_conventions={
        'shift_phi':-90,
        'flip_theta':True,
        'flip_phi':False,
        'in_degrees':True,
        'add_invisible_sky':True
    }

antenna_hpmap_inst = antenna_inst.convert2hp(frequency=45, **update_antenna_conventions)

# hp maps are by default in galactic coordinates, so we need to specify the LST and latitude of the local observer
lst = 18
LATITUDE = -35.206667
rotation_parameters = create_rotation_parameters(lst, LATITUDE)
rotator = create_rotator(lst, LATITUDE, coord=["G", "C"])

antenna_hpmap = antenna_hpmap_inst.get_map(rotator=rotator)

# galaxy
projview(
    antenna_hpmap,
    cmap='jet',
    return_only_data=False,
    graticule=True,
    graticule_labels=True,
    title='Galactic coordinates',
    xtick_label_color='w',
    # projection_type='cart'
)

# from galaxy to local
projview(
    antenna_hpmap,
    cmap='jet',
    return_only_data=False,
    coord=['G','C'],
    rot=rotation_parameters,
    graticule=True,
    graticule_labels=True,
    title='Local coordinates',
    xtick_label_color='w',
    # projection_type='cart'
)